% Natasha Savage, Maria Minakova


close all;
clear all;

tic

global k1a 
global k1b 
global k2a  
global k2b 
global k3 
global k4a 
global k4b 
global k5a 
global k5b 
global k6a 
global k6b 
global k7 
global eta
global dx
global dx2
global kbirth %Masha dev sept2012
global k8a %rate of pheromone binding to receptor complex                              %Masha dev sept2012
global k8b %rate of pheromone unbinding to receptor complex                            %Masha dev sept2012
global k9 %rate of nucleotide exchange in Cdc42D->T assisted by liganded receptor RecA %Masha dev sept2012

ic_conc     = [];    % holds the concentration in the internal compartment 
                     % [Cdc42Tic Cdc42Dic vSNARE,RecCmplxic] %Masha dev sept2012
cyt_Cdc42   = [];
 
Cdc42_store   = fopen('data_Cdc42_time_course', 'w'); 
neg_store     = fopen('data_negative',          'w'); % going to hold the time points and protein if it goes negative
TOP_store     = fopen('data_TOP_time_course',   'w'); % putting TOP info in her for plotting
mass_store    = fopen('data_mass',              'w'); % going to hold total mass
RecCmplx_store =fopen('data_receptor_time_course','w'); %storing time evolution of the receptor complex concentration on the membrane %Masha dev sept2012
RecA_store     =fopen('data_reca_time_course',  'w'); %storing time evolution of the receptor complex bound to pheromone on the membrane %Masha dev sept2012

% variables for time-stepping
dt=0.05;            % time step in seconds
maxt = 20/dt;       % number of time-steps in a 20 sec interval
%maxinterval = 60*60/20;   % number of 20sec intervals to simulate
maxinterval = 60*60/20;   % number of 20sec intervals to simulate %Masha debug

% variables for spatial discretization
N = 100;                            % number of spatial grid points
Cellsize = 5;                       % sphere diameter (microns)
Internalcompartment = 0.7;          % relative size (area) of internal compartment vs plasma membrane 
dx  = Cellsize*sqrt(pi)/N;          % spatial meshsize (1D)
dx2 = sqrt(Internalcompartment)*dx; % spatial meshsize for internal compartment (1D)
dx0 = dx;                           % record original dx
%Rc = 3; %micron according to the Gorychev's paper, consistent with eta0=0.01 %Masha dev sept2012

nt     = 0;     % counter
sec    = 0;     % counter for sampling
curr_t = 0;     % keeps track of current time

%% initalising diffusion mechanism ________________________________________

nRsteps = 100;      % take nRsteps reaction steps per diffusion step
dt2 = dt/nRsteps;   % reaction time-step

totcnt = 1;

%% reaction constants

mult  = 10;
mult2 = 2;

k1a = 10;       %s-1,       BemGEFc -> BemGEF
k1b = 10;       %s-1,       BemGEF -> BemGEFc

% GEF
k2a = 0.16;     %uM-1.s-1,  BemGEF + Cdc42D -> Cdc42T
k3  = 0.35;     %uM-1.s-1,  BemGEF42 + Cdc42D -> Cdc42T
%GAP
k2b = mult2*0.315;    %s-1,       Cdc42T -> Cdc42D 

k4a = 10;       %uM-1.s-1,  BemGEF + Cdc42T -> BemGEF42
k4b = 10;       %s-1,       BemGEF42 -> BemGEF + Cdc42T
k7  = 10;       %uM-1.s-1,  BemGEFc + Cdc42T -> BemGEF42

% Cdc42D on
k5a = mult*0.9;      %s-1,       GDI42c -> GDI42
k6b = mult*0.5;      %s-1,       GDI42 -> Cdc42D + GDIc
% Cdc42D off
k6a = mult*1.5;      %uM-1.s-1,  Cdc42D + GDIc -> GDI42
k5b = mult*0.13;     %s-1,       GDI42 -> GDI42c

%Receptor (Complex) birth and sink rates [microM/s] %Masha dev sept2012
kbirth = 0.0021143; % [uM/s] given 0.006*Vcell volume of yeast golgi and synthesis rate 4 molecules/s per cell. Cell radius is 5 micron.
%kbirth = 0.0 %no synthesis, for debugging: plot internalization with time, find halftime -> kbirth from there.
%RecCmplx + alpha <-> RecA %Masha dev sept2012
k8a=2.0 % binding of pheromone. Simon PNAS, 2003,100(19),10764-10769
k8b=0.01 %unbinding of pheromone. Simon PNAS, 2003,100(19),10764-10769

%Cdc42D + RecA ->(k9)-> Cdc42T + RecA
%k9=0.35; %nucleotide exchange facilitated by Receptor bound to pheromone complex
k9=0.0 %no effect on Cdc42, for debugging: plot internalization with time, find halftime -> kbirth from there.

eta0=0.01;  %eta: Vm/Vc, membrane/cytoplasm volume correction
eta =eta0;

% as size of outer membrane changes so does eta
mem_depth  = (Cellsize/2)*(1-(1/((eta0+1)^(1/3))));
cyt_mult   = ((Cellsize/2)^3)/(eta0+1);
Rnew_mult  = N/2/(pi^0.5);
eta = ((Rnew_mult^3 * dx^3) - (((Rnew_mult * dx)-mem_depth)^3))/cyt_mult;

% compute diffusion-coeff times dt/dx^2
Dconst = 0.0025;                % diffusion coefficient
Diff1  = Dconst * dt/dx^2;

% diffusion matrix, includes dt/dx^2 factor
e        = ones(N,1); 
Dxx      = spdiags([e, -2*e, e], -1:1, N, N);
Dxx(1,N) = 1;                   % periodic boundary conditions
Dxx(N,1) = 1;
Dxx      = Dxx*Diff1;

Mxx = kron(speye(N),Dxx);
Myy = kron(Dxx,speye(N));
Hop = speye(N^2) - Mxx - Myy;   % heat operator, I-Dxx

%% Initial conditions, homo steady state + perturbation

% initial concentrations
  Cdc42T                  = zeros(N); 
  Cdc42T(N/2    ,N/2    ) = 500; 
  Cdc42T(N/2    ,(N/2)+1) = 500; 
  Cdc42T((N/2)+1,N/2    ) = 500; 
  Cdc42T((N/2)+1,(N/2)+1) = 500; 
  BemGEF42                = zeros(N);
  BemGEF                  = zeros(N);
  BemGEFc                 = 0.017*ones(N);
  Cdc42D                  = zeros(N);
  GDI42                   = zeros(N);
  GDI42c                  = 5 - ((4*500)*eta/N/N)*ones(N);
  GDIc                    = ((4*500)*eta/N/N)*ones(N);

% initial trafficing concentrations, membrane and internal compartment
vSNARE   = 1.8015*ones(N);     % V-snare concentrations on membrane
vSNAREic = 1; % V-snare concentrations in the internal compartment (ic)
RecCmplx = 3.1778*ones(N); %Receptor initial concent on the memb. [microM], given cell rad 5micron, mem thickness 10 nm, tot number 6000 per cell %Masha dev sep2012
RecCmplxic = 0.0 %Receptor Complex initial concentration in the internal compartment (ic) %Masha dev sep2012
RecA = zeros(N); %Pheromone bound receptor, init conc [microM] on the membrane.

% homo steady state for 5uM system - doesn't matter as IC is accessible:
        % Cdc42T   = 15.0653
        % BemGEF42 = 1.1712
        % Cdc42D   = 22.5098
        % GDI42    = 29.5416
Cdc42Dic = 18.8832; 
GDI42ic  = 27.9147; 

GDI42c = GDI42c - eta*(GDI42ic + Cdc42Dic)*dx2*dx2/dx/dx;
GDIc   = GDIc   + eta*(          Cdc42Dic)*dx2*dx2/dx/dx;


%% initialize vesicle mechanism ___________________________________________

% Parameters related to vesicle fusion/fission

TOP = 100;              % number of bins in the window of highest Cdc42/GTP
                        % concentration
nVInlen  = N/50;        % no. of bins (1D) used to represent an exocytic 
                        % vesicle
nVOutlen = N/100;       % no. of bins (1D) used to represent an endocytic 
                        % vesicle
periodOut = 0.6;        % mean period for endocytosis, sec (= 15 sec actin 
                        % patch/25 patches in cell)
periodIn = 4*periodOut; % mean period for exocytosis, sec
patch_fill = 10;        % fill level for endocytosis, used for sink_case =2
patch_fill_give_up = 24;% even if the patch doesn't fill up with cargo, 
record_vSNARE = [];     % record [v-SNARE] at time of endocytosis
record_endot = [];      % record time of endocytosis
record_lifetime = [];   % record lifetime of each sink
%Masha dev sept2012
record_RecCmplx = [];     % record [RecCmplx] at time of endocytosis
record_RecA = [];     % record pheromone bound receptor [RecA] at time of endocytosis

endo_red_time = 15;     % 15 s before endocytosis when spots are "red"
exo_cnt = 0;            % count exocytic events
xg = (repmat(1:N,N,1)); yg = xg';  % xg, yg used to to identify window
           
% random seed
 RandStream.setGlobalStream(RandStream('mt19937ar','seed',0));
%rand('seed',0); %Masha debug RandStream does not exist in MAtlab 7.6
% Finding the initial TOP bins
Cdc42T_totm = Cdc42T+BemGEF42;                % total Cdc45/GTP membrane
                                              % bistribution

[val,Tbins] = sort(Cdc42T_totm(:),'descend'); % 'val' = sorted concentrations
                                              % 'Tbins' = sorted index    
circ_mem = Tbins(1:TOP)';


circ_ind    = [];               % will hold positions of TOP bins
ves_prob    = [];               % value 1 in TOP bins, 0 otherwise
ves_outprob = [];               % value 0 in TOP bins, 1 otherwise

Ratio = 40;         % probability of a bin endocytosing within the window
                    % are 40x higher than the probability of a endocytosing
                    % outside

% probability that endocytosis occurs inside TOP bins
endo_m = 40/139;

% storage for wrapping, cells chosen to endocytose, counters, and tracking Cdc42 distribution
Fwarp = zeros(N+1); % storage for use in interpolation that needs wrap around func values
endo_cell        = struct('endo_time',[],'pick_time',[],'vSNARE',[],'RecCmplx',[],'RecA',[],'Cdc42T',[], ...
    'Cdc42D',[], 'BemGEF42',[], 'GDI42',[], 'BemGEF',[], 'posx',[],'posy',[],'dx',[],'num',0); %added RecCmplx,RecA %Masha dev 2012
endo_cell_record = struct('endo_time',[],'pick_time',[],'vSNARE',[],'RecCmplx',[],'RecA',[],'Cdc42T',[], ...
    'Cdc42D',[], 'BemGEF42',[], 'GDI42',[], 'BemGEF',[], 'posx',[],'posy',[],'dx',[],'num',0); %added RecCmplx,RecA %Masha dev 2012

%% ________________________________________________________________________
% write initial data into files

      for i = 1:N
        fprintf(Cdc42_store,'%5.4f ',Cdc42T(i,:)+Cdc42D(i,:)+GDI42(i,:)+BemGEF42(i,:));
        fprintf(Cdc42_store,'\n');
        fprintf(RecA_store,'%5.4f ',RecA(i,:));
        fprintf(RecA_store,'\n');
        fprintf(RecCmplx_store,'%5.4f ',RecA(i,:));
        fprintf(RecCmplx_store,'\n');
      end
      
      cyt_Cdc42 = [cyt_Cdc42 GDI42c(1,1)];
      
      % Calclate the global mass __________________________________________
     Cdc42s = (mean(mean(GDI42c))) ...                      % Cdc42D cytoplasmic content 
     + eta*(mean(mean(GDI42+Cdc42D+Cdc42T+BemGEF42))) ...   % Cdc42 membrane content 
     + eta*(Cdc42Dic*(dx2^2)/(dx^2)) ...                    % Cdc42D on inner membrane 
     + eta*( GDI42ic*(dx2^2)/(dx^2)) ...                    % GDI42 on inner membrane           
     + eta*sum(endo_cell.Cdc42T(1:endo_cell.num)   .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % Cdc42T in vesicle                            
     + eta*sum(endo_cell.Cdc42D(1:endo_cell.num)   .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % Cdc42D in vesicle 
     + eta*sum(endo_cell.GDI42(1:endo_cell.num)    .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % GDI42 in vesicle
     + eta*sum(endo_cell.BemGEF42(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);        % BemGEF42 in vesicle 

     Bem1s = (mean(mean(BemGEFc))) ...                      % cytoplasmic content 
     + eta*(mean(mean(BemGEF+BemGEF42))) ...                % Bem1 membrane content  
     + eta*sum(endo_cell.BemGEF(1:endo_cell.num)   .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % Bem1 vesicle content 
     + eta*sum(endo_cell.BemGEF42(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);        % BemGEF42 vesicle content                                                        

     GDIs = (mean(mean(GDIc+GDI42c))) ...                   % GDI cytoplasmic content 
     + eta*( GDI42ic*(dx2^2)/(dx^2)) ...                    % GDI42 on inner membrane 
     + eta*(mean(mean(GDI42))) ...                          % GDI on outer membrane
     + eta*sum(endo_cell.GDI42(1:endo_cell.num).* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);            % GDI42 in vesicle

     vSNAREs        = eta*mean(mean(vSNARE)) ...              % v-SNARE membrane content
     + eta*(vSNAREic*(dx2^2)/(dx^2)) ...                   % v-SNARE on inner membrane content 
     + eta*sum(endo_cell.vSNARE(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);          % v-SNARE vesicle content 
 
     RecCmplxs        = eta*mean(mean(RecCmplx)) ...              % Receptor complex membrane content %Masha dev sept2012
     + eta*(RecCmplxic*(dx2^2)/(dx^2)) ...                   % Receptor complex on inner membrane content 
     + eta*sum(endo_cell.RecCmplx(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);          % Receptor complex vesicle content 
     
     RecAs        = eta*mean(mean(RecA)) ...              % Receptor complex membrane content %Masha dev sept2012 
     + eta*sum(endo_cell.RecA(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);          % Receptor complex vesicle content 

      %____________________________________________________________________                                                  
     
      % writing global mass into file _____________________________________
      fprintf(mass_store,'%5.4f\t%5.4f\t%5.4f\t%5.4f\t%5.4f\t%5.4f\n',Cdc42s, Bem1s, GDIs,vSNAREs,RecCmplxs,RecAs); %added Receptore complex storing %Masha dev sept2012

      % calculate the internal compartment mass ___________________________
      % ic_conc = [Cdc42Tic GDI42ic vSNARE,RecCmplxic] %Masha dev sept2012
      ic_conc(totcnt,1) = (Cdc42Dic);   % Cdc42D  on inner membrane 
      ic_conc(totcnt,2) = (GDI42ic);
      ic_conc(totcnt,3) = (vSNAREic);   % v-SNARE on inner membrane         
      ic_conc(totcnt,4) = (RecCmplxic); %Receptor complex on inner membrane %Masha dev sept 2012
      %____________________________________________________________________
            
      totcnt = totcnt+1;
%%
%Masha debug
frcount=0;

for intervals=1:maxinterval     % each loop is 20 seconds
  
  for t0=1:1:maxt               % dt = 0.05 incraments: maxt = 20/dt
    
    t = (intervals-1)*maxt+t0;
    curr_t = curr_t+dt;         % update clock
    
    react_step_Euler_VSS_GDI42ic_tolAdjust
    
    
    % Finding the TOP bins after each react_step
    Cdc42T_totm = Cdc42T+BemGEF42;                % total Cdc45/GTP membrane
                                                  % bistribution
                                            
    [val,Tbins] = sort(Cdc42T_totm(:),'descend'); % 'val' = sorted concentrations
                                                  % 'Tbins' = sorted index    
    circ_mem = Tbins(1:TOP)';
                
    circ_ind = circ_mem;                          % index of TOP concentrations

    ves_prob = zeros(N); ves_prob(circ_ind)=1;    % 1 in TOP bins, 0 otherwise
    ves_outprob = ones(N);                        % 0 in TOP bins, 1 otherwise
    ves_outprob(circ_ind)=0;

    % plotting to check this is working OK - - - - - - - - - - - - - - - - 
    % writing plot info into a file for checking at end of simulation
    [x,y] = ind2sub([N,N],circ_ind);
    if(t>=1/dt && mod(t,1/dt)==0)   % store every second
        save output.mat
        for i = 1:N
            fprintf(TOP_store,'%5.4f ',Cdc42T_totm(i,:));
            fprintf(TOP_store,'\n');
        end
        fprintf(TOP_store,'%5.4f ',x);
        fprintf(TOP_store,'\n');
        fprintf(TOP_store,'%5.4f ',y);
        fprintf(TOP_store,'\n');
        fprintf(TOP_store,'%5.4f ',Cdc42T_totm(circ_ind));
        fprintf(TOP_store,'\n');
%         fprintf(RecCmplx_store,'%3.2f ',t*dt); %Masha dev sept2012
%         fprintf(RecCmplx_store,'%3.3f',mean(mean(RecCmplx))); %Masha dev sept2012
%         fprintf(RecCmplx_store,'\n'); %Masha dev sept2012
%         fprintf(RecA_store,'%3.2f ',t*dt); %Masha dev sept2012
%         fprintf(RecA_store,'%3.3f',mean(mean(RecA))); %Masha dev sept2012
%         fprintf(RecA_store,'\n'); %Masha dev sept2012
    end
    % END plotting to check this is working OK - - - - - - - - - - - - - -                 

    cargo_step_WT_GDI42ic       

    %Masha debug
%     if(t>=3/dt && mod(t,3/dt)==0) % take a mov frame 10 seconds _________
%      frcount=frcount+1;
%      %fig1=figure(1);
%      %image(Cdc42T+GDI42+BemGEF42);
%      xlim([0 N]);
%      ylim([0 N]);
%      title('GTP bound forms of Cdc42: Cdc42T + GDI42 + BemGEF42');
%      %cdc42fr(frcount)=getframe(fig1);
%      %cdc42fr= addframe(cdc42fr,fig1);
%      fr=getframe(fig1);
%      writeVideo(cdc42fr,fr);
%      fig2=figure(2);
%      image(RecA);
%      title('Receptor bound to pheromone');
%      %recafr(frcount)=getframe(fig2);
%      %recafr= addframe(recafr,fig2);
%      fr=getframe(fig2);
%      writeVideo(recafr,fr);
%      fig3=figure(3);
%      image(RecCmplx);
%      title('Receptor complex not bound to pheromone');
%      %reccfr(frcount)=getframe(fig3);
%      %reccfr= addframe(reccfr,fig3);
%      fr=getframe(fig3);
%      writeVideo(reccfr,fr);
%     end
    
    %if(t>=1/dt && mod(t,1/dt)==0) % take a snap shot every second _________
    if(t>=5/dt && mod(t,5/dt)==0) % take a snap shot every 5 seconds _________ %Masha debug
      save output.mat
      for i = 1:N
        fprintf(Cdc42_store,'%5.4f ',Cdc42T(i,:)+Cdc42D(i,:)+GDI42(i,:)+BemGEF42(i,:));
        fprintf(Cdc42_store,'\n');
        fprintf(RecA_store,'%5.4f ',RecA(i,:));
        fprintf(RecA_store,'\n');
        fprintf(RecCmplx_store,'%5.4f ',RecCmplx(i,:));
        fprintf(RecCmplx_store,'\n');
      end
      
      cyt_Cdc42 = [cyt_Cdc42 GDI42c(1,1)];
      
      % Calclate the global mass __________________________________________
     Cdc42s = (mean(mean(GDI42c))) ...                      % Cdc42D cytoplasmic content 
     + eta*(mean(mean(GDI42+Cdc42D+Cdc42T+BemGEF42))) ...   % Cdc42 membrane content 
     + eta*(Cdc42Dic*(dx2^2)/(dx^2)) ...                    % Cdc42D on inner membrane 
     + eta*( GDI42ic*(dx2^2)/(dx^2)) ...                    % GDI42 on inner membrane           
     + eta*sum(endo_cell.Cdc42T(1:endo_cell.num)   .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % Cdc42T in vesicle                            
     + eta*sum(endo_cell.Cdc42D(1:endo_cell.num)   .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % Cdc42D in vesicle 
     + eta*sum(endo_cell.GDI42(1:endo_cell.num)    .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % GDI42 in vesicle
     + eta*sum(endo_cell.BemGEF42(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);        % BemGEF42 in vesicle 

     Bem1s = (mean(mean(BemGEFc))) ...                      % cytoplasmic content 
     + eta*(mean(mean(BemGEF+BemGEF42))) ...                % Bem1 membrane content  
     + eta*sum(endo_cell.BemGEF(1:endo_cell.num)   .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % Bem1 vesicle content 
     + eta*sum(endo_cell.BemGEF42(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);        % BemGEF42 vesicle content                                                        

     GDIs = (mean(mean(GDIc+GDI42c))) ...                   % GDI cytoplasmic content 
     + eta*( GDI42ic*(dx2^2)/(dx^2)) ...                    % GDI42 on inner membrane 
     + eta*(mean(mean(GDI42))) ...                          % GDI on outer membrane
     + eta*sum(endo_cell.GDI42(1:endo_cell.num).* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);            % GDI42 in vesicle

     vSNAREs        = eta*mean(mean(vSNARE)) ...              % v-SNARE membrane content
     + eta*(vSNAREic*(dx2^2)/(dx^2)) ...                   % v-SNARE on inner membrane content 
     + eta*sum(endo_cell.vSNARE(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);          % v-SNARE vesicle content 
     
     RecCmplxs        = eta*mean(mean(RecCmplx)) ...              % Receptor complex membrane content %Masha dev sept2012
     + eta*(RecCmplxic*(dx2^2)/(dx^2)) ...                   % Receptor complex on inner membrane content 
     + eta*sum(endo_cell.RecCmplx(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);          % Receptor complex vesicle content 
    
     RecAs        = eta*mean(mean(RecA)) ...              % Receptor complex membrane content %Masha dev sept2012
     + eta*sum(endo_cell.RecA(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);          % Receptor complex vesicle content 
     %____________________________________________________________________                                                  
     
      % writing global mass into file _____________________________________
      fprintf(mass_store,'%5.4f\t%5.4f\t%5.4f\t%5.4f\t%5.4f\n',Cdc42s, Bem1s, GDIs,vSNAREs,RecCmplxs,RecAs);

      % calculate the internal compartment mass ___________________________
      % ic_conc = [Cdc42Tic GDI42ic vSNARE]
      ic_conc(totcnt,1) = (Cdc42Dic);   % Cdc42D  on inner membrane 
      ic_conc(totcnt,2) = (GDI42ic);
      ic_conc(totcnt,3) = (vSNAREic);   % v-SNARE on inner membrane                                      
      ic_conc(totcnt,4) = (RecCmplxic);  %Receptor complex on inner membrane %Masha dev sept 2012     
      %____________________________________________________________________
            
      totcnt = totcnt+1;
      RecCmplxic; %Masha debug
      %mean(mean(RecA)) %Masha debug
    end
          
  end  

end

%% ________________________________________________________________________

fclose(Cdc42_store);
fclose(neg_store);
fclose(TOP_store);
fclose(mass_store);
fclose(RecCmplx_store); %Masha dev sept 2012
save vSNAREic_end.mat vSNAREic
save RecCmplxic_end.mat RecCmplxic %Save last snapshot of Receptor complex on inner membrane %Masha dev sept 2012
save GDI42ic_end.mat GDI42ic
save Cdc42Dic_end.mat Cdc42Dic

save vSNARE_end.mat vSNARE
save RecCmplx_end.mat RecCmplx %Save last snapshot of Receptor complex on outer membrane %Masha dev sept 2012
save RecA_end.mat RecA %Save last snapshot of Receptor bound to pheromone on outer membrane %Masha dev sept 2012
save GDI42_end.mat GDI42
save Cdc42D_end.mat Cdc42D
save Cdc42T_end.mat Cdc42T
save BemGEF_end.mat BemGEF
save BemGEF42_end.mat BemGEF42

save GDIc_end.mat GDIc
save GDI42c_end.mat GDI42c
save BemGEFc_end.mat BemGEFc

time_whole = toc/60/60
save output.mat

